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ABSTRACT 

Chemistry has a key role in the evolution of the interstellar medium (ISM), so it 
is highly desirable to follow its evolution in numerical simulations. However, it may 
easily dominate the computational cost when applied to large systems. In this paper 
we discuss two approaches to reduce these costs: (i) based on computational strategies, 
and (ii) based on the properties and on the topology of the chemical network. The 
first methods are more robust, while the second are meant to be giving important 
information on the structure of large, complex networks. 

To this aim we first discuss the numerical solvers for integrating the system of ordinary 
differential equations (ODE) associated with the chemical network. We then propose 
a buffer method that decreases the computational time spent in solving the ODE 
system. We further discuss a flux-based method that allows one to determine and 
then cut on the fly the less active reactions. 

In addition we also present a topological approach for selecting the most probable 
species that will be active during the chemical evolution, thus gaining information on 
the chemical network that otherwise would be difficult to retrieve. This topological 
technique can also be used as an a priori reduction method for any size network. 
We implemented these methods into a ID Lagrangian hydrodynamical code to test 
their effects: both classes lead to large computational speed-ups, ranging from x2 
to x5. We have also tested some hybrid approaches finding that coupling the fiux 
method with a buffer strategy gives the best trade-off between robustness and speed- 
up of calculations. 
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fragmentation of the interstellar gas (|jappsen et al.l [20051 ) 
and it is one of the most powerful diagnostic way to relate 
with the observational signatures coming from the last 
generation of telescopes like ALMAQ and Herschefl Unfor- 
tunately tracking its evolution needs large computational 
resources. The complexity of the problem arises from its 
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mathematical representation: a network of chemical reac- 
tions has its counterpart in a system of coupled ordinary 
differential equations (ODE) that is often stiff (i.e. with 
coefficient highly- varying in size). The solution of such 
a system may becomes too expensive computationally, 
especially when the number of reactions is large, or when 
the system is formed by many elements of fluid. 

Over the years several approaches have been pro- 
posed to solve this problem: reducing the number of reac- 
tions/species making physical or chemical assumptions is 
one of the most used. The number and the kind of reactions 
being removed depends on the as trophysical environmen t 
that one wants to simulate (e.g. iNelson fc Langeil Il999l ). 
This method can be useful when the network is small, but it 
is hard to apply when a larger network is employed, leading 
to large errors with complex networks. 

Another widely used appro ach is to reduce th e number 
of non-equilibrium species (e.g. iGlover et al.]|2010l ) treating 



2 T. Grassi et al. 



the remaining species as being in equilibrium. This method 
decreases the chemical timestep and thus the computational 
cost, but it is problem-dependent since the selection of the 
non-equilibrium species depends on the environment one 
wants to simulate. 

In this paper we discuss a set of methods that can be 
used to increase the feasibility of accurate simulations based 
on large chemical networks while also providing some key 
information on the structure of the network. 

We introduce the problem in Sect [21 then we consider 
two classes of reduction strategies: (i) the first is based on 
methods more related to computer science and we call them 
"pure" computational strategies (see Sect[3|, (ii) while the 
second set is based on the analysis of the properties and on 
the topology of the chemical network (SectQ. In Sect [5] we 
will show tests on some of the methods previously discussed. 



2 OVERVIEW OF THE PROBLEM 

A chemical network is composed of a number of reactions 
each generally represented byA-|-B — > C-l-D, and the prob- 
ability that a reaction occurs is related to a rate coefficient 
ki{T) that depends on the gas temperature T. In this frame- 
work each species evolves following an ordinary differential 
equation (ODE) that takes into account the reactions that 
form and destroy the ith species 

^ kj{T)nriQ)nr2U)-ni ^ kj{T)nri(j) , (1) 

jGform j^destr 

where Mri(j) and ?^r2(j) are the number densities of the first 
and the second reactants of the jih reaction, and the sub- 
script form and destr represents the set of the reactions 
that form and destroy the ith species. Note that for three- 
body reactions we must multiply the first part (and/or the 
second, depending on the reaction) of the RHS term by the 
number density of the third reactants, i.e. w^sjj). The net- 
work is then described by a system of ODEs. 

From a computational point of view the problem of solv- 
ing a system of ODEs resides mainly in building the RHS 
term of Eqn.Q, while the total number of equations has a 
smaller (but n on-negligible) impact on the solver efficiency 
(|Tuppei|[2003 '). To achieve this we can either remove the 
species, decreasing the number of ODEs and hence the di- 
mensionality of the problem, and/or reduce the number of 
reactions, thus diminishing the time spent to build the RHS 
of Eqn.([T|. Both approaches will be discussed in Sect 2] 



3 PURE COMPUTATIONAL STRATEGIES 

In this Section we discuss the strategies that are less related 
to the properties and the features of the chemical network, 
but allow to improve the computational performance of solv- 
ing a system of ODEs aimed at tracking the chemical evolu- 
tion. We first present (i) an analysis of the different solvers, 
then (ii) the buffering method, and finally (iii) a general 
discussion on standard computational techniques. 



3.1 Numerical solvers 

Astrochemical networks are represented by a system of stiff 
ODEs, which requires appropriate solvers. Due to the nu- 
merical instability of the system, one of the most used solvers 
is the implicit backward differencing (BDF) scheme, but 
depending on its implementation it can greatly affect the 
computational performance. A large number of solvers have 
been proposed over the years like Rosenbrock, Bulirsch- 
Stoe r-Bader-Deufl hard, implicit Runge-Kutta, and Gears 
(Pre ss et al.lll99^ ). but unfortunately these methods have 
poor performance when applied to astrochemical networks 
since a large number of iterations and function evalua- 
tions are required to integrate the ODE system. A good 
way to solve this problem is to use the well-established 
ODEPACK/SUNDIALS packageS (|Hindmars h et al. 200^), 
a solver suite based on GEARS and LSODE that allows 
to solve stiff ODE systems more efficiently. The DVODE 
fits our class of problems as well, but nevertheless, when 
the system present s a sparse Jacobian the DLSODES solver 
(|Hindmarshlll983l '). which has the capability of handling 
sparse matrices, i s more effici ent. DLSODES has been al- 
ready discussed in lNeiadl(|2005l ). and we have tested hero the 
two solvers on a larger network jWakelam fc He rbst 200^) 
involving 452 species and more than 4000 reactions. The 
DVODE and the DLSODES will be tested in SectlST] 



3.2 Buffering 

In the framework of hydrodynamical simulations one usually 
solves the chemical ODEs for each fluid element (i.e. each 
gas particle or grid cell), the number of which can easily 
exceed 10®, depending on the resolution: for this reason re- 
ducing the number of calls to the solver saves a large amount 
of CPU-time. This allows us to take advantage of the fact 
that we can avoid to solve the chemical ODEs for fluid el- 
ements which have chemical conditions (e.g. temperature, 
species numerical densities, . . . ) similar to those that were 
already calculated, since they will show the same evolution 
with time. This comes from the fact that a system of ODEs 
is a Cauchy's problem that depends on the set of initial con- 
ditions and on the differential equations which are weighted 
by the values of the reaction rates. 

The buffering method we propose consists of storing the 
already calculated chemical results. If the initial conditions 
of a given fluid element matches the initial conditions of a 
buffered element we can avoid to call the solver again, and 
use instead the stored evolution. In the opposite case we call 
the solver to evolve the fluid element and we store the initial 
and the flnal conditions in the buffer (see also Algorithm [T] 
discussed below). 

To determine if two particles are similar we loop over 
the initial conditions and we check if the expression Si = 

—Xi\/xi < e is true, where Xij is the ith initial condition 
of the jth particle of the buffer, and Xi is the ith initial 
condition of the test particle. It is worth noting that we 
also need a constraint on the time-steps of the two particles, 
namely |df — dt^l < et, where et is the time threshold, and in 
the tests of Sect l5.2l we use et ~ 0.1 yr. The time threshold is 

^ computation.llnl.gov/casc /sundials / main.html 



Reduction of astrochemical networks 3 



Algorithm 1 - Pseudocode for the buffering method proposed here. Note that the term particle is the same as fluid element. 
The aim of this psuedocode is to find the final conditions x of a particle with initial conditions x. See details in Sect l3.2l 



1 


X initial conditions 


// X are the initial conditions of a given particle 


2 


for (j = A^, 1) do 


/ / loop over butter particles {j — N ^ 1) 


3 


found-f- True 


/ / set default value of found ttag 


4 


if {\dtj - dt\ > et) then 


/ / compare dt 


5 


found-<— False 


II -i--C J-O J- J?! 

/ / set found ttag to false 


6 


skip to next fluid element 


/ / skip to the next burter particle (next j) 


7 


end if 


1 1 


8 


for (i — 1, Nc) do 


1 1 ^ • 'J. • 1 J 'J. • ( ■\ 

1 1 loop over initial conditions \i) 


9 


if [\xi —Xji\/Xi > e) then 


/ / check smnlarity tor the ith initial condition 


10 


found-<— False 


// if in J ri 

/ / set found ttag to false 


11 


exit from initial condition loop 


111 1 1 ■■j*i i^i* /'\ 

/ / breaks loop over initial conditions [%) 


12 


end if 


1 1 

II 


13 


end for 


1 1 end loop over initial conditions 


14 


if (found) then 


1 1 particles match 


15 


X ^ Xj 


/ / load final conditions from buffer 


16 


exit from buffer loop 


// go to the next cell (break loop j) 


17 


end if 


II 


18 


end for 


/ / end loop over buffer particles 


19 


if (not found) then 


//if the particle is not in the buffer 


20 


X -s— solver(x, dt) 


/ / compute final conditions with solver 


21 


XiV+l = X 


/ / store final conditions in the buffer 


22 


Xjv+1 = X 


/ / store inital conditions in the buffer 


23 


dtN+1 = dt 


/ / store time step in the buffer 


24 


N ^ N + 1 


/ / increase the size of the buffer 


25 


end if 


// 



chosen in order to have a good approximation over the time- 
step and for this reason its value must be a fraction of the 
mean hydrodynamical time-step, depending on the desired 
accuracy. From our test we found that a small et provides 
accurate results even with larger values of e. 

The total error is given by the chosen e and et, and 
it also depends on the chaotic behaviour of the chemical 
network. If the system is unstable to small perturbations of 
the initial parameters we must choose a very small e thereby 
reducing the advantage of using the buffer. 

In order to increase the speed-up of the method we scan 
the buffer from the most recent values to those stored at 
earlier stages, since it is more likely that the last particles 
stored are similar to the test particle. This increases the 
probability of finding a similar particle earlier during the 
scan, allowing us to break the cycle after a few comparisons: 
this increases the advantage of using this method. 

Since the size A'' of the buffer increases while the sim- 
ulation is running, we must define a maximum buffer size, 
namely Nmax, that is A''max = A'tot • m where A'tot is the 
total number of fluid elements employed in the simulation, 
and m ^ 1 (we set m = 20 and A^tot ~ fOO in our tests). 
When the buffer is full (A'^ = A'max) we simply move the last 
A'max/S buffered cells at the beginning of the buffer and we 
set A^ = A'max /2, in order to re-use the last A''max /2 cells 
stored. 

We provide a pseudocode (Algorithm[T]) in order to clar- 
ify the above description. Note that the pseudocode refers to 
a given fluid element (called test element) with initial con- 
ditions X (line 1) that is compared to all the particles in the 
buffer via the loop j = TV — s- 1. The aim of this algorithm is 
to find the final conditions x for the test particle by compar- 



ison, (line 15) or, if comparison fails, via a call to the solver 
(line 20). In the algorithm we first check the condition for 
the time-step dt (line 4) that is written in order to be false 
when the time-steps of the test particle and the jth buffer 
particle are different (line 6) to avoid useless calculations 
(i.e. skip to the next buffer particle). Then it loops over the 
set of particles to compare the initial conditions of the test 
particle and of the jth buffer particle (line 8). In this case 
we also write the conditions in order to save CPU time by 
exiting from the initial conditions loop (line 11). If both of 
these tests (time-step and initial conditions) are false, the 
variable found remains set to true (line 14), and in this case 
the jth buffer particle is similar to the test particle: we load 
then the final conditions from the buffer particle (line 15) 
and we break the loop over the buffer (line 16). If, after the 
loop over the buffer, the variable found (line 19) is still false 
we need to calculate the final conditions in the standard way 
(via solver, line 20), we store the new particle in the buffer 
(lines 21 and 22), and we increase the size of the buffer (line 
24). 

In the tests proposed in Sect[Sl 70% of the calls to the 
solver are saved with a speed-up of x 5 and a very small error 
on the chemical evolution. We also found, as expected, that 
the number of skipped calls to the solver decreases during 
the simulation, because when the shock is fully developed 
there are less similar fluid elements than at the beginning, 
but, more in general, the total number of evolved fluid ele- 
ments loaded from the buffer (and hence not calculated with 
the solver) depends on the features of the simulated envi- 
ronment. Note that a large buffer (Amax S> A'ceii) reduces 
the rate of efficiency decreasing, altough a large buffer needs 
a non- negligible computational time to be scanned. 
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3.3 Other strategies 

There are other strategies related to standard programming 
techniques that can give an additional speed-up as are the 
various optimization flags allowed by widely-used compil- 
ers. We don't discuss here the details which are compiler- 
dependent, but we remind the reader that aggressive opti- 
mization (often indicated as -03 flag), inter-procedural op- 
timization, and loop unroUing/vectorization can lead to sig- 
nificant speed-ups. We also remark that good programming 
practices (e.g. save divisions, avoid casting, pre-calculate 
mostly used expressions, . . . ) help the global performances 
of the code and should be always employed. These strategies 
coupled with an accurate code profiling will make the code 
considerably less CPU-demanding. 

Evaluating the rate coefficients is another source of 
CPU-consumption. Often the ODE for the temperature is 
coupled with the others ODEs through cooling and/or heat- 
ing and trough the adiabatic index (7) that depends on the 
mass fractions of the species. For this reason the rate co- 
efficients must be evaluated at each solver's step, and usu- 
ally these coefficients contain intrinsic functions as exp() and 
log() that are expensive in terms of computer time. One way 
to solve this problem (while keeping the same accuracy) is 
to tabulate the rates and to interpolate them during the 
simulation when required, obtaining a speed-up up to a x5 
depending on the reaction rates used. 

The last solver-related issue concerns the function called 
by the solver that contains the differential equations sur- 
mised by Eqn.(IT|. These equations can be either explicitly 
written equation-by-equation in the code, or called by using 
a loop. The latter approach produces a more compact code 
which is easy to handle and modify, and it is faster during 
compilation (especially for very large networks), but unfor- 
tunately it is less efficient at runtime. In the tests proposed 
here we have preferred for the sake of simplicity the loop ap- 
proach, even if explicitly writing the equations would give 
an additional speed-up. 



4 THE REDUCTION METHODS 

In this Section we review those reduction strategies that are 
more directly related to the properties and the features of 
the network itself. We briefiy discuss the most used meth- 
ods for reducing the workload associated with the chemistry, 
and we also introduce a new approach based on the topology 
of the network. The need to reduce the astrochemical net- 
works arises from the huge computational time required for 
astrophysical simulations. Common practices, based on ar- 
bitrary cutting, are often used to deal with this problem, e.g. 
those based on chemical or physical considerations. As al- 
ready stated in the Introduction, this approach is reasonable 
for very small networks where the complexity of the connec- 
tions among the species is clear enough to decide what to 
remove. 

There are several methods aiming at reducing the com- 
putational cost of the chemistry. A common practice is to 
pre-calculate the models with some one- zone chemical code 
in order to create a database of model evolutions, and then, 
when needed, to interpolate the final conditions during a 
large simulation. This method works fine for a small num- 



Table 1. Overview of the methods discussed in this Paper, where 
d.o.f. means degrees of freedom, transf. transformation, and infos 
is for informations on the structure of the chemical network. See 
text for further details. 



Method 


pros 


cons 


Direct 


exact 


very slow 


Interp. 


fast 


few d.o.f. 


ANN 


fast, many d.o.f. 


hard fine-tuning 


User-based 


fast 


arbitrary 


Lumping 


fast 


small-systems 


SVD, PCA 


fast 


transf. overhead 


Buffer 


fast 


efficiency decreases 


Flux-based 


fast, infos 


DLSODES overhead 


Topology 


fast, infos 


misses initial cond. 



ber of parameters, e.g. the temperature and the initial num- 
ber densities of the species involved. When the number of 
species grows, the dimensionality of the database increases 
thus making the database huge and difficult to interpolate. 

The latter problem can be solved with an artificial neu- 
ral network (ANN) approach which works similarly to an 
interp olator, but can han dle a larger number of free param- 
eters jCrassi et aLllioill ). Unfortunately an ANN needs a 
lot of fine-tuning to reproduce the results of the database, 
and for this reason most of the times it cannot be included 
in large simulations as an interpolator. 

Another class of reduction m ethods is based on linear 
algebra, (i) The lumping method (|Okino fc Mavrovounioti^ 
11998) allows one to reduce the dimensionality of the problem 
by making some linear combination of the species involved. 
This method has been developed for small systems and for 
this reason cannot be applied to complex networks, (ii) Also 
the singular value decomposition (SVD) and the principal 
component analysis (PCA) are aimed at reducing the dimen- 
sionality of the ODE system using a transformation matrix 
that allows to redefine the coordinates of the space of the 
parameters (in our case temperature and all the numerical 
abundances) in order to define a new space where the di- 
mensionality of the problem is smaller. The main drawback 
is that computing this transformation matrix has a non- 
negligible computational cost and it could easily result in a 
computational overhead. 

The overview of these methods is presented in Tab[TJ 
which includes their pros and the cons features in terms of 
usage. 

We tackle this computational problem by employing in- 
stead a mor e robust techn ique known as the fiux -method 
presented in lTupped l|2002l ) and lGrassi et al.l (|2012l '). It is an 
on the fly procedure aimed at determining the less active re- 
actions and then reducing the n umber of terms i n th e RHS of 
Eqn.(IT|. The tests proposed in lCrassi et al.l (|2012l ) showed 
large speed-ups ranging from x2 to xlO depending on the 
network employed. It is worth noting that this method has 
been developed for the solver DVODE, while when applied 
to DLSODES solver it generates an overhead since it in- 
terferes with the solver's subroutines aimed at handling the 
sparsity of the Jacobian (Tupper, private communication). 

A new a priori method based on the connectivity fea- 
tures of a network (i.e. on topology) is therefore presented 
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in the next Section, and some of its key aspects and results 
will be discussed there. 

4.1 The topology of astrochemical networks 

An astrochemical network can be viewed as a directed graph 
where the vertexes (or nodes) are the species of the chem- 
ical system, while the edges (or links) are determined by 
reactions between species. A reaction as A + B — ^ C is rep- 
resented with three nodes (namely A, B, C) and two edges 
(A ^ C and B — >■ C). When we consider astrophysical net- 
works the number of species/vertexes spans from a few tens 
to more than several hundreds, while reactions/edges can 
be more than five thousands. This large number of items 
becomes complicated to be displayed as a graph and, as a 
consequence, the topological features are usually not evident 
at first sight. Different methods have been proposed through 
the years to explore the properties of the n etworks as sum- 
marized bv iJollev fc Douglas! (|2010l . |2012| ). A widely used 
method consists in measuring the degree of the nodes that 
comprise the network, where degree is defined as the sum of 
the connection of a node to its neighbourhoods. This quan- 
tity shows some interesting properties when applied to real 
examples. One of the most intriguing is the power-law dis- 
tribution of the probability of finding a node with a given 
number of connections. More explicitly the probability of 
finding a node with a degree d decreases as P{d) oc 
(where 7 > is a parameter), it means that finding a node 
with many connections is less probable that finding a less 
connected one. 

This property is true for scale-free networks, that are 
cons iderably robust and stable to random node removing 
(e.g. lBarabaslll999t ). In particular when the network is ro- 
bust it allows to remove some nodes without damaging the 
global stability of the network because the largest part of 
the information fiows through the so called hubs that are 
the most connected elements of the system. When we deal 
with astrophysical networks we are interested in finding sta- 
ble structures from a "noisy background" represented by the 
less active nodes or, as in the scheme just depicted, we want 
to know what node can be removed or not. 

Even if astrochemical networks are not scale-free, since 
they have an exponential degree distribution instead of a 
power law l|jollev fc Douglasll2010l ). we can determine which 
are the sub-structures that can be deleted by using some 
ranking techniques. The degree of a node can be a good way 
to determine whether a node is important or not, and then 
if it could be removed or not; when we decide to use the 
degree as a criterion we assume that it is very probable that 
a large part of the information will flow through very con- 
nected nodes instead of flowing through the less connected. 
It is evident that removing vertexes with many connections 
(hubs) could seriously damage the whole network because of 
their crucial role in the global network activity. 

Another approach to rank the most important species is 
to calculate how much each node is linked to nodes that are 
highly connected to the rest of the network. This method 
results from the assumption that the most connected nodes 
carry the largest part of the global network information and 
in particular we choose that being the neighbour of a very 
connected node increases the probability of becoming active 
during the network evolution. This fact can become more 



clear when we take a social network as an example: in this 
case a person (a node) can easily reach critical information 
if it is connected to a few individuals that massively partic- 
ipate to the information exchange {hubs), rather than being 
in contact with a large number of poorly connected persons. 

We call this "second degree" and it can be represented 
as the sum of the degrees of the nodes close to a given node 
(i.e. its neighbours). For the ith node is 

(^)d. = ^ , (2) 

where '■^^d; is the second degree of the ith node, Ni is the 
set of the neighbours of the ith node, and dj is the degree 
of the node j. 

In this way the second degree becomes a parameter that 
is more powerful than the degree itself, because it takes into 
account not only the activity of a given node, but also the 
activity of its close neighbours. In principle we can also cal- 
culate the nth degree by iteratively applying the idea above 
for n ^ 00 

<"^d. = . (3) 

This is similar to the Google PageRank (IPage et al.lll999h 
that assigns a numerical weighting to each web page de- 
termining its relative importance. Both algorithms have an 
iterative approach aimed at finding the probability distri- 
bution that represents the likelihood that a given species 
participates to the global network activity, or a person ran- 
domly clicking on various links will arrive at a particular 
page. 

In this first analysis, to compute the second degree we 
consider the graph as undirected, since to determine the ac- 
tivity of a given node we have to include the edges that leave 
that node, which represent the destruction reactions, and 
the edges pointing at that node, being reactions which form 
the species represented by the given node. Under this as- 
sumption the activity of the node is determined by both the 
reactions that form and destroy the corresponding species, 
hence the need to consider the graph as undirected. 

There are three approaches to calculate the second de- 
gree and the difference between them is determined by the 
weight assumed during the calculation of the degree di in 
the equation 

where Wj is the weight. The first method to obtain di is sim- 
ply (i) to count the number of connections that reach or leave 
a node assuming Wj — 1. This method allows to compute 
the second degree without considering the time evolution, 
but it neglects whether the reactions carrying information 
to the node are important or not. (ii) To improve on that 
choice, one can change the weight to Wj = kij where kij 
is the rate coefficient of the reaction that involves the ith 
node as product and the jth node as a reactant. This is 
more accurate than the first approach, but it cannot be cal- 
culated once and for all since kij is temperature-dependent. 
This method also ignores the fact that at a given time a 
reaction could not be active because the abundance of the 
reactants is zero, (iii) The latter issue can be avoided by 
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considering Wj = Fir where Fij is the flux of a reaction cal- 
culated as Fij = kij rij rir where rij is the abundance of 
the neighbour reacting species, and the product runs over 
the other reactants excluding rij. This calculation can be 
only performed at runtime, because one needs to know the 
abundances of the species during the simulation, and more- 
over, the rate coefflcient kij depends on the temperature, 
which changes during the evolution of the model. 

In this paper we shall show the results only for the first 
method that can be used as a good estimator to determine 
the most important species, allowing to eliminate from the 
network the less active ones before calculation. The other 
two methods have more general validity and accuracy but re- 
quire evaluation during runtime, which is a significant draw- 
back from the computational standpoint. 

Once we have ranked the species in the network we can 
choose a criterion to remove the less important ones, and 
to choose this we want to make sure that we do not re- 
move species with high abundances. To cope with this risk 
we introduce a second degree threshold d„ based on species 
densities 

'■^'dn — min ['^'dij Vi | > nt . (5) 

that depends on the lowest of the second degrees ^^^di of the 
species with a number density rii greater than a user-defined 
density threshold nj. We also define a second threshold df to 
keep important species when the first threshold can be ex- 
cessively constraining (e.g. initial monoatomic gas), namely 

^""^df = /-max[''>d,] . (6) 

which is determined by a certain fraction / ^ 1 of the max- 
imum second degree. 

We decide to cut the species that do not satisfy both 
criteria at the same time, hence the final threshold becomes 

(^)d, = minp)d„,(^)d,] , (7) 

which gives the final criteria '^'d^ > ^-^^dt employed to decide 
whether or not the ith species is kept in our calculation. In 
the tests we present in Sect l5.3l we use / = 0.6 and nt — 
10"^ 



the tests presented in this paper employ the DLSODES 
solver with an absolute and relative tolerance of 10"'"^ and 
lO"^'^ respectively (which are accurate enough compared 
to the values of the number densities usually employed in 
ISM simulation), and an inter nally-generated Ja cobian (op- 
tion MF=222, more details in lHindmarshlll983l ). Note also 
that the code has largely been optimized and compiled with 
Intel® Fortran Composer XE 2011 Update 6 using a stan- 
dard set of optimization flags. 

It is important to remark that since the scope of 
these tests is to determine the computational efficiency of 
the methods proposed, we use the OSU database as it is 
(osu_01_2007, 4431 reactions and 452 species) and, more- 
over, we do not include any cooling nor heating. Under the 
latter assumption the hydrodynamics of the shock is not 
affected by the chemical evolution. However, the tests dis- 
cussed here represent - computationally speaking - a typical 
ISM scenario with a large chemical network. We plan to de- 
scribe the coupling between hydrodynamics, chemistry and 
cooling together with reduction methods in a forthcoming 
paper. 

The methods are all compared with the full evolution, 
i.e. without any reduction. To show the accuracy of the dif- 
ferent methods we have plotted the abundances of a given 
species in all the gas shells for each reduction method com- 
pared to the same shell in the full model. More in detail 
the plots represent the accuracy in reproducing the chemi- 
cal behaviour of the gas during the evolution of the shock in 
all the shells for a given species. If a method reproduces ex- 
actly the chemical evolution of the full model all the points 
will lie along a straight line (i.e. y = x), while the distance 
from the line increases if the method fails to reproduce the 
original model. These plots represents the dispersion of the 
values obtained for a given method (rimothod) compared to 
the values of the full model nfun. We have also included two 
lines representing the magnitude of fiuctuations of one order 
of magnitude in each direction (i.e. y = 10 a; and y = 0.1 a; 
labelled -|-odm and -odm respectively). 

We found that coupling together buffer and flux meth- 
ods gives the best results both in terms of robustness and of 
speed-up. 



5 TEST CASES 

We propose here tests involving a selection of the aforemen- 
tioned metho ds, empl o ying the OSU reactions database as 
in|Wakclam fc HerbstI (|2008l ') - without PAHs - coupled to a 
simple ID Lagrangian code aimed at simulating a Sedov-like 
gas shock (jBodenheimer et alllioo^ . The gas is composed 
of a set of spherical shells representing an inner region with 
initial conditions Tin = lO'' K and pin = lO"'^'^ g/crr? , and 
another set of shells for the outer region with Tout = 10 K 
and Pout = 10~^^ g/cm^ and radius R = Ipc. The number 
of the inner shells is Nin = 40, while the number of the 
outer shells is = 60, giving Mot = 100. Note the results 
found for the tests presented in this Section are still valid 
for larger systems and also for systems based on an Eulerian 
approach, and 3D hydrodynamics. 

The initial abundances a re the same as the EA2 model 
of IWakelam fc HerbstI l|2008l ) (i.e. solar metallicity initial 
conditions) and we let the system evolve for 10*' yr. All 



5.1 Comparison of solvers 

The characteristics of the chemical network play a key role 
to determine what is the best numerical solver. In this frame- 
work we propose a ID shock model test mainly focused 
on Jacobian spars i ty, using DVODE and D LSODES solvers 
(|Hindmarshlll983l : iHindmarsh et al.l 120051 ). It is important 
to remark that employing an internally generated or a user- 
provided Jacobian will affect the performance of the solver, 
since building the Jacobian has a non- negligible cost for the 
solver; however for the sake of simplicity in all our tests we 
use this latter method. 

The Jacobian associated with the network of the OSU 
database is extremely sparse (~ 94%), which is very com- 
mon for astrochemical networks. For this reason in the shock 
test framework the numerical solver DLSODES outperforms 
the more general DVODE solver achieving approximately a 
xlOO speed-up. 
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5.2 Buffer method 

We test the buffer method with two thresholds, namely 
e — 10^^ and e — 10~^, the latter being a finer approxima- 
tion. We found good results especially for the 10~® threshold 
which reproduces the full model almost exactly, but, as ex- 
pected, when increasing the threshold value (e — 10~^) the 
method starts to be less accurate in reproducing the origi- 
nal data, but the error remains always below one order of 
magnitude (see Fig[T]and Fig[2}. 



5.3 Topological method 

The topological method reproduces the full model with a 
good accuracy even if some dispersion is observed for the 
and H species (Fig[T] top and FiglU bottom). The threshold 
we chose selects 2369 reactions over 4431, hence it removes 
almost 50% of the network connections. As stated in the pre- 
vious Section, the reduction based on the topological method 
has a probabilistic meaning and its best property lies in the 
capability of giving information on which species/reactions 
are important or not. For instance, our tests show that the 
most important hubs (the most connected species) are, in 
that order, free electrons, H, H2, He^, C, H"*", O and C""", as 
expected. This topological approach finds that species like 
CI, Mg, Fe, and their ions, but also MgH, and HF and many 
C-chains species, are not so important within this network, 
thus many reactions involving He^ and C"*" are neglected. 
The latter can explain the fluctuation of species as shown 
in Fig[T] We want to underline once again that, unlike flux 
or buffer, the topological method is mainly aimed at giving 
a priori chemical information on the network rather than 
strongly reducing the number of the reactions, or the species. 
A more robust reduction technique based on topology has al- 
ready been discussed: using the flux as the weight in Eqn.(|4]) 
will provide a more accurate analysis. However, this method 
is not a priori, and for this reason it could generate a large 
overhead when the reduction is less effective. Moreover, from 
a computational point of view, this technique is closer to the 
flux method, since the weights are exactly the fluxes and the 
two methods are largely comparable. The only advantage is 
the amount of topological information provided during the 
system evolution. Due to these considerations, it seems likely 
that flux-based reduction techniques are preferable when a 
robust reduction method is required, even if the error asso- 
ciated with the topology method is generally less than one 
order of magnitude. Depending on the desired accuracy, the 
method can nevertheless be useful in some applications. 



5.4 The flux-based method 

For the flux method presented here we assume that the 
length of the hydrodynamical time-step is equal to the 
length of the macro-step, which is a good approximation 
since temperature and total density will remain constant 
during a hydrodynamical time-step. Under this assumption 
the evaluation of the fluxes is made at the beginning of each 
hydrodynamical time-step for each gas shell. This allows one 
to couple the DLSODES solver with the flux method. The 
results (Fig[T]and Fig[2} are in good agreement with the full 
solution and can be further improved using a more accurate 



Table 2. Normalized CPU-time measured for different reduction 
methods. See text for further details. 



Method <CPU 



Full 1.00 

Topology 0.37 

Buffer (-1) 0.40 

Buffer (-5) 0.61 

Flux 0.50 

Flux-I-Buffer 0.24 



Flux-I-Buffer-I-Topology 0.17 



Table 3. Speed-up overview for the methods discussed here. Note 
that the values listed here refer to the ISM model presented in 
this paper. One can find different speed-ups depending on the 
chemical network employed, the numerical framework, and some 
other parameters. However, these values are representative for the 
methods discussed and can be used as reference. 



Method 


Speed 


-up 


Solvers 


xlOO 




Reduction methods 


up to 


xlO 


Buffering 


up to 


x5 


Rates tabulation 


x5 





definition of macro-step (see iGrassi et al.l I2OI2I for further 
details) . 

5.5 Hybrid methods 

The last method discussed here is a mixing between buffer 
(e = 10~^), flux, and topology. The accuracy of this hybrid 
method is determined by the least accurate approach as we 
can see in Fig[l](top), where topology fails, and also in Fig[2] 
(bottom) where topology is not as accurate as the other 
methods. 

The normalized CPU-times of the different methods are 
reported in Table (2] as expected hybrid methods are the 
fastest, but coupling flux and buffer shows a better accuracy 
than the coupling of flux, buffer, and topology. We also ob- 
tain a large speed-up for the coarsest of the two buffer meth- 
ods, while the one with e = 10~^ has the smallest speed-up, 
although it has the best accuracy among the methods pro- 
posed. Finally, the topology shows good speed-up with less 
accurate results. It is worth noticing that the latter reduc- 
tion method gives important information about the chemical 
network due to its topological approach, and that to really 
see the effects of such a cut on the global hydrodynamical 
evolutions we need a test which couples the chemistry to the 
dynamics via heating and cooling, and via photon diffusion. 



6 CONCLUSIONS 

In this Paper we have discussed some techniques aimed at 
reducing the computational cost of including a chemical net- 
work into astrophysical codes that simulate the evolution of 
the ISM. We have divided the methods into two classes, 
namely the "pure" computational strategies that are more 
related to computer science, and the reduction methods that 
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are based on the properties and on the structure of the chem- 
ical network. In the first class we discussed the solvers em- 
ployed, the buffering method, and a brief review on other 
strategies, while in the second class we described the flux 
method and introduced a topology-based method. 

To test the different methods we employed a ID hy- 
drodynamical Sedov-like shock test with a standard OSU 
(osu_0I_2007) chemical network, which is a large chemical 
network (> 4000 reactions). The methods tested are buffer, 
topology, and flux, including some hybrid methods that cou- 
ple the previous ones. These methods reproduce the re- 
sults of the full model with good speed-ups, except for the 
topology-based approach which has a larger error (compared 
to the other methods), since it is mainly devoted to rank by 
importance the different species rather than providing a ro- 
bust reduction technique as the flux method. Nevertheless, 
the results provided by this method suggest that the topo- 
logical approach can lead to an efficient reduction based on 
the features of the chemical network, which gives at the same 
time some critical information about the global properties 
of the interconnections between the various species included 
in the ISM model. 

Finally, we found that the most robust method among 
those examined turns out to be the coupling between the flux 
approach and the buffer technique, and this hybrid method 
also achieves one of the best speed-ups (almost x5) as sug- 
gested by our tests. The fastest shock simulation is obtained, 
as expected, by using flux, buffer, and topology all together, 
but the topology affects the global behaviour of this hybrid 
method producing results that are less accurate than flux 
and buffer methods and their hybrid. 

A final overview of the speed-ups obtained from the 
various methods presented in this paper is given in Table [S] 
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Figure 1. Comparison of the full model with the different re- 
duction methods for 0+ (top) and CO (bottom). As indicated in 
the key the methods are topology (dots), buffer with e = 10~^ 
(triangles up) and e = 10~^ (open circles), flux (triangles down), 
and hybrid (diamonds). For colour figure see the online version. 
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Figure 2. Comparison of the full model with the different re- 
duction methods for OH (top) and H (bottom). As indicated in 
the key the methods are topology (dots), buffer with e = 10~^ 
(triangles up) and e = 10~^ (open circles), flux (triangles down), 
and hybrid (diamonds). For colour figure see the online version. 
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